function [STAT,tb] = pupil_permutation(c1, c2, nRep)
%
% Takes two conditions, c1 and c2, and finds time regions of significant
% difference. c1 and c2 should be cell arrays of individual participant
% averages traces of the same length.
%
% NB: The user is entrusted with ensuring that all traces (for each
% participant in each condition) are the same length. There are no checks
% for this, but if it is not the case, unexpected results can arise.
%
% "nRep" is the number of iterations in the permuation test.

fs = 500;
nSs = length(c1);
if nSs ~= length(c2)
    error('Different numbers of subjects in the two conditions')
end

datalen = length(c1{1}); % assumes all lengths equal!

%% Organize data into matrices for simpler calculation
c1 = cell2mat(c1');
c2 = cell2mat(c2');

%% T-Tests at each time point to get T-values from actual data
actual_diff = zeros(1, datalen);
for t = 1:datalen
    actual_diff(t) = fast_Ttest(c1(:,t), c2(:,t));
end

%% Permutation test: Iteratively re-label, re-calculate T-values
c1_shuf = c1;
c2_shuf = c2;
all_perms = zeros(nRep, datalen);
disp('starting permutation test')
for p = 1:nRep
    % Randomly re-label
    [c1_shuf, c2_shuf] = mix_conditions(c1_shuf, c2_shuf);
    % Re-calculate the T-Test
    for t = 1:datalen
        all_perms(p,t) = fast_Ttest(c1_shuf(:,t), c2_shuf(:,t));
    end  
end

%% Where is the T-value from the real data > from the permuted at p < 0.05?
STAT = zeros(1,datalen);
for t = 1:datalen
    STAT(t) = sum(all_perms(:,t) >= actual_diff(t)) / nRep;
end

tb = (1:datalen) ./ fs;

end


%% --- Condition shuffler

function [x_shuf, y_shuf] = mix_conditions(x, y)

which2shuf = logical(datasample([0,1], size(x, 1)));
x_shuf = x; y_shuf = y;
x_shuf(which2shuf,:) = y(which2shuf,:);
y_shuf(which2shuf,:) = x(which2shuf,:);

end

%% --- Faster t-statistic?

function T = fast_Ttest(x, y)

T = abs(mean(x - y) / (std(x - y)/sqrt(length(x))));

end










